1. upload MPA locations
  2. link MPA locations to grid cells on globe
MPAs <- read.csv(paste0(data_dir, "notakeonly_coord.csv"))
head(MPAs)
##   Centroid_Longitude Centroid_Latitude
## 1           177.8167         -17.35000
## 2           -17.0500          20.78333
## 3           179.1040         -17.10960
## 4           179.9667         -16.85000
## 5          -159.7500         -18.86670
## 6           158.1257           6.79908
tail(MPAs)
##     Centroid_Longitude Centroid_Latitude
## 304           8.559806          42.35926
## 305           3.164926          42.47173
## 306           6.390493          43.00589
## 307         -81.396041          19.39510
## 308        -111.091496          25.85519
## 309         -42.750000         -62.75000
plot(MPAs$Centroid_Longitude, MPAs$Centroid_Latitude)

# Transform longitude from MPAs
MPAs$Centroid_Longitude_trans <- MPAs$Centroid_Longitude
MPAs$Centroid_Longitude_trans[which(MPAs$Centroid_Longitude<0)] <- 360 + MPAs$Centroid_Longitude[which(MPAs$Centroid_Longitude<0)]


dfresults <- read.csv(paste0(results_dir,"SigmaD_Hem.csv"))
dfresults <- dfresults[order(dfresults$No),]
head(dfresults)
##       No NN.sigma_2100_4.5_today NN.station_2100_4.5_today
## 17414  1            0.0008453082                     30792
## 17415  2            0.0003516396                     30792
## 17416  3            0.0018296424                     30792
## 17417  4            0.0032007127                     30792
## 17418  5            0.0074244659                     30575
## 17419  6            0.0055364983                     30575
##       NN.Mdist_2100_4.5_today numPCs_2100_4.5_today NN.sigma_today_2100_4.5
## 17414              0.03673379                     2                2.745675
## 17415              0.02368996                     2                2.692024
## 17416              0.05405385                     2                2.677252
## 17417              0.07151309                     2                2.705799
## 17418              0.10900851                     2                2.776901
## 17419              0.09409835                     2                2.891746
##       NN.station_today_2100_4.5 NN.Mdist_today_2100_4.5 numPCs_today_2100_4.5
## 17414                     40066                3.196740                     2
## 17415                     40066                3.145595                     2
## 17416                     40066                3.131513                     2
## 17417                     40066                3.158725                     2
## 17418                     40066                3.226513                     2
## 17419                     40066                3.336049                     2
##       NN.sigma_2100_8.5_today NN.station_2100_8.5_today NN.Mdist_2100_8.5_today
## 17414                5.866236                     30680                6.201382
## 17415                5.993259                     30680                6.324892
## 17416                6.079694                     30680                6.408981
## 17417                6.122091                     30680                6.450239
## 17418                6.114727                     30680                6.443073
## 17419                6.062258                     30680                6.392015
##       numPCs_2100_8.5_today NN.sigma_today_2100_8.5 NN.station_today_2100_8.5
## 17414                     2                8.160708                     40066
## 17415                     2                8.292361                     40066
## 17416                     2                8.160708                     40066
## 17417                     2                8.209536                     40066
## 17418                     2                8.209536                     40066
## 17419                     2                8.292361                     40066
##       NN.Mdist_today_2100_8.5 numPCs_today_2100_8.5   lat long
## 17414                8.512412                     2 -69.5 20.5
## 17415                8.470758                     2 -68.5 20.5
## 17416                8.465380                     2 -67.5 20.5
## 17417                8.501935                     2 -66.5 20.5
## 17418                8.577215                     2 -65.5 20.5
## 17419                8.693851                     2 -64.5 20.5
##       NN.sigma_1800_today NN.station_1800_today NN.Mdist_1800_today
## 17414         0.005340623                 40066           0.0924152
## 17415         0.007501424                 40009           0.1095737
## 17416         0.009629160                 40067           0.1241975
## 17417         0.015490226                 40009           0.1577081
## 17418         0.024201364                 40066           0.1974682
## 17419         0.049860481                 40066           0.2848815
##       numPCs_1800_today NN.sigma_today_1800 NN.station_today_1800
## 17414                 2         0.019575390                 26505
## 17415                 2         0.020981196                 29981
## 17416                 2         0.014089519                 26082
## 17417                 2         0.006465863                 26082
## 17418                 2         0.011645698                 25749
## 17419                 2         0.005112661                 25973
##       NN.Mdist_today_1800 numPCs_today_1800 lat_1800_today long_1800_today
## 17414          0.17743251                 2          -69.5           379.5
## 17415          0.18374455                 2          -69.5           378.5
## 17416          0.15036689                 2          -68.5           379.5
## 17417          0.10170870                 2          -69.5           378.5
## 17418          0.13663930                 2          -69.5           379.5
## 17419          0.09041724                 2          -69.5           379.5
##       latshift_1800_today latshift lat_today_2100_8.5 long_today_2100_8.5
## 17414                  NA        0              -69.5               379.5
## 17415                  NA       -1              -69.5               379.5
## 17416                  NA       -1              -69.5               379.5
## 17417                  NA       -3              -69.5               379.5
## 17418                  NA       -4              -69.5               379.5
## 17419                  NA       -5              -69.5               379.5
##       lat_today_2100_4.5 long_today_2100_4.5
## 17414              -69.5               379.5
## 17415              -69.5               379.5
## 17416              -69.5               379.5
## 17417              -69.5               379.5
## 17418              -69.5               379.5
## 17419              -69.5               379.5
tail(dfresults)
##          No NN.sigma_2100_4.5_today NN.station_2100_4.5_today
## 17408 40115                4.501010                       554
## 17409 40116                4.525272                       554
## 17410 40117                4.569456                       554
## 17411 40118                4.606448                       554
## 17412 40119                4.625980                       554
## 17413 40120                4.632551                       554
##       NN.Mdist_2100_4.5_today numPCs_2100_4.5_today NN.sigma_today_2100_4.5
## 17408                4.879349                     2            5.414412e-05
## 17409                4.902748                     2            7.803501e-06
## 17410                4.945368                     2            5.824597e-04
## 17411                4.981062                     2            3.656808e-04
## 17412                4.999911                     2            3.734047e-03
## 17413                5.006253                     2            1.033971e-02
##       NN.station_today_2100_4.5 NN.Mdist_today_2100_4.5 numPCs_today_2100_4.5
## 17408                      8267             0.009295342                     2
## 17409                     24925             0.003528828                     2
## 17410                     12303             0.030490772                     2
## 17411                     21911             0.024158378                     2
## 17412                     16363             0.077249962                     2
## 17413                     16363             0.128716445                     2
##       NN.sigma_2100_8.5_today NN.station_2100_8.5_today NN.Mdist_2100_8.5_today
## 17408                8.292361                     40112                9.377801
## 17409                8.292361                     40113                9.465370
## 17410                8.292361                     40113                9.475441
## 17411                8.292361                     40113                9.455885
## 17412                8.292361                     40113                9.425748
## 17413                8.292361                     40113                9.404554
##       numPCs_2100_8.5_today NN.sigma_today_2100_8.5 NN.station_today_2100_8.5
## 17408                     2               0.9093877                      9460
## 17409                     2               1.2192752                      8629
## 17410                     2               1.4372294                      8629
## 17411                     2               1.6185109                      8629
## 17412                     2               1.7344557                      8629
## 17413                     2               1.7884077                      8629
##       NN.Mdist_today_2100_8.5 numPCs_today_2100_8.5  lat  long
## 17408                1.423342                     2 84.5 379.5
## 17409                1.733062                     2 85.5 379.5
## 17410                1.945650                     2 86.5 379.5
## 17411                2.120635                     2 87.5 379.5
## 17412                2.231984                     2 88.5 379.5
## 17413                2.283687                     2 89.5 379.5
##       NN.sigma_1800_today NN.station_1800_today NN.Mdist_1800_today
## 17408        0.0033782225                 37019          0.07347198
## 17409        0.0041649559                 37769          0.08159263
## 17410        0.0001037297                  3229          0.01286606
## 17411        0.0047175381                  7551          0.08684629
## 17412        0.0377153031                  7551          0.24717272
## 17413        0.0646634655                  7551          0.32537576
##       numPCs_1800_today NN.sigma_today_1800 NN.station_today_1800
## 17408                 2            1.728258                   613
## 17409                 2            1.931682                   613
## 17410                 2            2.170640                   613
## 17411                 2            2.350300                   554
## 17412                 2            2.424026                  3223
## 17413                 2            2.455741                  3336
##       NN.Mdist_today_1800 numPCs_today_1800 lat_1800_today long_1800_today
## 17408            2.226041                 2           86.5           347.5
## 17409            2.420731                 2           86.5           353.5
## 17410            2.648753                 2           88.5            59.5
## 17411            2.819972                 2           88.5           103.5
## 17412            2.890216                 2           88.5           103.5
## 17413            2.920432                 2           88.5           103.5
##       latshift_1800_today latshift lat_today_2100_8.5 long_today_2100_8.5
## 17408                  -2       -2               71.5               130.5
## 17409                  -1       -1               73.5               119.5
## 17410                  -2       -2               73.5               119.5
## 17411                  -1       -1               73.5               119.5
## 17412                   0        0               73.5               119.5
## 17413                   1        1               73.5               119.5
##       lat_today_2100_4.5 long_today_2100_4.5
## 17408               77.5               113.5
## 17409               83.5               239.5
## 17410               81.5               157.5
## 17411               79.5               218.5
## 17412               78.5               183.5
## 17413               78.5               183.5
world <- map_data("world2")
 ggplot()+
   geom_polygon(data = world, aes(x=long, y = lat, group=group)) + 
   geom_point(data = MPAs, aes(x= (Centroid_Longitude_trans), y = Centroid_Latitude)) + 
   stat_summary_2d(data=dfresults, aes(x=long, y = lat, z= NN.sigma_today_2100_4.5), bins=80, alpha = 0.8)

# Transform longitude from this analysis
dfresults$long_trans <- dfresults$long
dfresults$long_trans[which(dfresults$long>360)] <- dfresults$long[which(dfresults$long>360)] - 360
  # 380 should correspond to 20; 360 corresponds to 0

 ggplot()+
   # geom_polygon(data = world, aes(x=long, y = lat, group=group)) + 
   geom_point(data = MPAs, aes(x= (Centroid_Longitude_trans), y = Centroid_Latitude)) + 
   stat_summary_2d(data=dfresults, aes(x=long_trans, y = lat, z= NN.sigma_today_2100_4.5), bins=80, alpha = 0.8)

summary(MPAs$Centroid_Longitude) # - 180 to 180
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## -177.37 -121.95  -76.01  -10.41  143.40  179.97
summary(dfresults$long) # 20 to 380
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    20.5   138.5   207.5   208.8   295.5   379.5
dfresults$MPA_nnIndex <- NA
dfresults$MPA <- FALSE
dfresults$MPA_nnDist <- NA
MPAs$No <- NA

for (i in 1:nrow(MPAs)){
  q <- as.matrix(cbind(MPAs$Centroid_Latitude[i], MPAs$Centroid_Longitude_trans[i])) # lat, long
  
  d <- as.matrix(dfresults[,c("lat", "long_trans")])
  nn <- get.knnx(d, q, k=1)
  
  MPAs$No[i] <- nn$nn.index
  MPAs$nnDist[i] <- nn$nn.dist
  MPAs$No_lat[i] <- dfresults$lat[nn$nn.index]
  MPAs$No_long_trans[i] <- dfresults$long_trans[nn$nn.index]
  #q
  #dfresults[nn$nn.index,c("lat", "long")]
  
  if(nn$nn.dist < 0.9){ # only include if station is close
    dfresults$MPA_nnDist[nn$nn.index] <- nn$nn.dist
    dfresults$MPA_nnIndex[nn$nn.index] <- nn$nn.index
    dfresults$MPA[nn$nn.index] <- TRUE
    }
}

# total number of MPAs included
length(which(MPAs$nnDist<0.9))
## [1] 177
# number of stations representing MPAs
sum(!is.na(dfresults$MPA_nnIndex))
## [1] 94
table(MPAs$No[which(MPAs$nnDist<0.9)])
## 
##   591  1098  1447  2657  4396  4856  8019  8091  8393  8452  8718  9503  9838 
##     1     1     1     3     1     1     1     1     1     1     3     1     1 
##  9839 10004 10088 10367 10369 10454 10530 10549 10550 10719 10825 10828 10830 
##     3     1     1     1     1     1     3     1     3     2     1     1     2 
## 10940 10946 11057 11059 11172 11180 11181 11188 11196 11439 11440 11441 11442 
##     1     1     1     1     1     2     2     1     1     1     1     2     3 
## 12168 12169 12272 12274 12534 12741 13594 13803 14377 15029 15157 15188 15189 
##     5     2     1     2     2     1     1     1     1     3     1     1     1 
## 15320 15350 15351 15885 16678 16833 17038 17197 18342 18765 18768 18920 18962 
##     1     3     1     1     1     3     1     1     1     1     3     1     2 
## 19095 19114 19266 19267 19571 23996 24133 24135 24263 24264 24390 24513 24634 
##     2     5     2     2     1     1     1     1     6    10     7     4     3 
## 24635 26464 27761 28671 28842 29090 29238 29512 29546 29665 29730 29731 29875 
##     7     2     1     2     1     1     1     1     8     2     3     1     1 
## 32070 33374 35810 
##     1     1     1
  #This should sum to 177
sum(table(MPAs$No[which(MPAs$nnDist<0.9)]))
## [1] 177
hist(dfresults$MPA_nnDist)

  1. look at degree of novelty for MPAs relative to non-MPA sites
ggplot(dfresults, aes(NN.sigma_today_2100_4.5,
                      fill = MPA)) + 
  geom_density(color="#e9ecef", alpha=0.3, position = 'identity') + theme_classic()

ggplot(dfresults, aes(NN.sigma_today_2100_8.5,
                      fill = MPA)) + 
  geom_density(color="#e9ecef", alpha=0.3, position = 'identity') + theme_classic()

  1. make arrow plots from 1800 to today to 2100 for MPAs

To Do 2: 1. Create fake environment envelope 2. Create fake locations with high ICV and different degrees of novelty 3. Create fake locations with low ICV and different degrees of novelty 4. Visualize

Visualize arrows

There is a confusing issue with this MPA analysis. I can’t just take the MPA column from dfresults and have it make sense all the time.

the 1800_today analysis takes a site today (given by lat,long) and compares it to the 1800 baseline. The NN in the 1800 baseline is lat_1800_today, long_1800_today. So if I want to draw an arrow from where a climate was in 1800 to where it is today today, I would draw from lat_1800_today to lat

the today_2100 analysis takes a site in 2100 (given by lat, long) and compares it to a today baseline. The NN in the today baseline is lat_today_2100, long_today_2100. So if I want to draw an arrow from where a climate is today in today’s baseline to where it will be in 2100, I would draw from lat_today_2100 to lat

ggplot() +
    geom_polygon(data = world, aes(x=long, y = lat, group=group)) + 
  # Plot from 1800 to location today
  # Arrows point from 1800 climate to MPA today
  geom_segment(data = dfresults[dfresults$MPA,],
               aes(x = long_1800_today, 
                   y = lat_1800_today, 
                    xend = long, yend = lat, 
                    colour = NN.sigma_1800_today), 
                 arrow = arrow(length = unit(0.1,"cm")))

ggplot() +
    geom_polygon(data = world, aes(x=long, y = lat, group=group)) + 
  # Plot from today to 2100 4.5
  # In 2100, what is the nearest neighbor today?
  # Arrow points from today's environment (long_today_2100_4.5) to the MPA in the future (long)
    geom_segment(data = dfresults[dfresults$MPA,],
               aes(xend = long, 
                   yend = lat, 
                  x = long_today_2100_4.5, 
                   y = lat_today_2100_4.5, 
                    colour = NN.sigma_today_2100_4.5), 
                 arrow = arrow(length = unit(0.2,"cm"))) 

ggplot() +
    geom_polygon(data = world, aes(x=long, y = lat, group=group)) + 
  # Plot from today to 2100 8.5
    geom_segment(data = dfresults[dfresults$MPA,],
               aes(xend = long, 
                   yend = lat, 
                  x = long_today_2100_8.5, 
                   y = lat_today_2100_8.5, 
                    colour = NN.sigma_today_2100_8.5), 
                 arrow = arrow(length = unit(0.2,"cm"))) 

Compare latitudes

poly <- data.frame(y=c(-90,90,90, -90, -90), x=c(-90,90,0,0,-90))


# If we consider a location in 2000 (lat), where is it's nearest neighbor 1800? (lat_1800_today)?
ggplot(dfresults) + geom_polygon(data=poly,aes(x,y), fill=adjustcolor("orange", 0.3)) +
   geom_point(aes(y=lat, x=lat_1800_today,
                                    color=NN.sigma_1800_today), alpha=0.2) +  
  scale_color_gradient2(low="red", high="blue", mid="grey", limits=c(0,8)) +
  theme_classic() + ylab("Latitude in 2000") +
  xlab("Latitude of nearest neighbor in 1800") + geom_abline(intercept=0,slope=1) 

# 2000 to 2100
# If we consider a location in 2100, where was the climate it came from in 2000?
ggplot(dfresults) +  geom_polygon(data=poly,aes(x,y), fill=adjustcolor("orange", 0.3)) +
  geom_point(aes(y=lat, x=lat_today_2100_4.5,
                 color=NN.sigma_today_2100_4.5), 
             alpha=0.2) +
  scale_color_gradient2(low="red", 
                        high="blue", mid="grey", limits=c(0,8)) +
  theme_classic() + ylab("Latitude in 2100 RCP 4.5") +
  xlab("Latitude of nearest neighbor in 2000") + geom_abline(intercept=0,slope=1) 

ggplot(dfresults) + geom_polygon(data=poly,aes(x,y), fill=adjustcolor("orange", 0.3)) +
  geom_point(aes(y=lat, x=lat_today_2100_8.5,
                 color=NN.sigma_today_2100_8.5), 
             alpha=0.2) +
  scale_color_gradient2(low="red", 
                        high="blue", mid="grey", limits=c(0,8)) +
  theme_classic() + ylab("Latitude in 2100 RCP 8.5") +
  xlab("Latitude of nearest neighbor in 2000") + geom_abline(intercept=0,slope=1) 

Latitudinal shifts

## Future latitude - past latitude
## 1800 to today
## "lat" is latitude of query in 2000
## lat_1800 is latitude of NN in 1800
dfresults$latshift_1800_today <- NA
condN <- which(dfresults$lat>0)
dfresults$latshift_1800_today[condN] <- dfresults$lat[condN] - dfresults$lat_1800_today[condN]
  # If 50 degrees in 1800 moved to 90 degrees in 2100 this
  # would give 90-50 = 40
  
condS <- which(dfresults$lat<=0)
dfresults$latshift_1800_today[condS] <- -1*(dfresults$lat[condS]-dfresults$lat_1800_today[condS])
  # If -50 degrees in 1800 moved to -90 degrees in 2100 this
  # would give -1*(-90--50) = 40


## Today to 2100 RCP 4.5
## Future latitude - past latitude
dfresults$latshift_today_2100_4.5 <- NA
dfresults$latshift_today_2100_4.5[condN] <- dfresults$lat[condN] - dfresults$lat_today_2100_4.5[condN]  
dfresults$latshift_today_2100_4.5[condS] <- -1*(dfresults$lat[condS] - dfresults$lat_today_2100_4.5[condS])

## Today to 2100 RCP 8.5
dfresults$latshift_today_2100_8.5 <- NA
dfresults$latshift_today_2100_8.5[condN] <- dfresults$lat[condN] - dfresults$lat_today_2100_8.5[condN]
dfresults$latshift_today_2100_8.5[condS] <- -1*( dfresults$lat[condS] - dfresults$lat_today_2100_8.5[condS])


par(mar=c(4,4,1,1), mfrow=c(3,1))
hist(dfresults$latshift_1800_today, col=adjustcolor("blue", 0.5), main="", xlab="Poleward shift from\nneareast neighbor in the past", breaks=seq(-90, 90, 1), xlim=c(-90,90))

hist(dfresults$latshift_today_2100_4.5, 
     col=adjustcolor("red", 0.2), main="", xlab="", 
     breaks=seq(-90, 90, 1), xlim=c(-90,90))

hist(dfresults$latshift_today_2100_8.5, 
     col=adjustcolor("green", 0.2), main="", xlab="",, 
     breaks=seq(-90, 90, 1), xlim=c(-90,90))

# color scale: if use "lat" it is the latitude in the later time
# color scale: if use "lat_1800_today" it is the latitude in the earlier time
# the x-axis is the shift of the latitude to that location
ggplot(dfresults) + geom_point(aes(x=latshift_1800_today, y = NN.sigma_1800_today, colour = lat), alpha=0.1)

ggplot(dfresults[dfresults$lat_1800_today>85,]) + geom_point(aes(x=latshift_1800_today, y = NN.sigma_1800_today, colour = lat), alpha=0.1)

# high latitudes in 1800 can only stay in the same place or move southward - color is latitude in 2000

ggplot(dfresults[dfresults$lat_1800_today>70 & dfresults$lat_1800_today<85,]) + geom_point(aes(x=latshift_1800_today, y = NN.sigma_1800_today, colour = lat), alpha=0.1)

ggplot(dfresults[dfresults$lat_1800_today<60 & dfresults$lat_1800_today>40,]) + geom_point(aes(x=latshift_1800_today, y = NN.sigma_1800_today, colour = lat_1800_today), alpha=0.1)

# 40-60 degress latitude in 1800 more northward


ggplot(dfresults) + geom_point(aes(x=latshift_today_2100_4.5, y = NN.sigma_today_2100_4.5, colour = lat_today_2100_4.5), alpha=0.1)

# colored by latitude in 2000, shows which way it shifted

ggplot(dfresults) + geom_point(aes(x=latshift_today_2100_8.5, y = NN.sigma_today_2100_8.5, colour = lat_today_2100_4.5), alpha=0.1)

# The RCP 8.5 shows that there aren't a lot of poleward shifts - because climates are so novel, the NN calculation becomes meaningless

# I wonder if some of the really negative latitudinal shifts are driven by distortion of the PC axes caused by the data not being multivariate normal.